# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/6/18

library(gplots)
library(Hmisc)
library(optparse)
library(ropls)
library(magrittr)
library(ComplexHeatmap)
library(tools)
library(tidyverse)

option_list <- list(
make_option("--i", default = "AllMet.csv", type = "character", help = "metabolite data file"),
make_option("--g", default = "true_SampleInfo.csv", type = "character", help = "sample group file"),
make_option("--et", default = "", type = "character", help = "extra data dir"),
make_option("--confounder", default = "", type = "character", help = "confounder data file"),
make_option("--config", default = "config.csv", type = "character", help = "config file")
)
opt <- parse_args(OptionParser(option_list = option_list))

createWhenNoExist <- function(f){
    ! dir.exists(f) && dir.create(f)
}

zd_pcor <- function (left = left , right =right  , confounding = confounding , unique_id_colname = "ID" , mannual_n_for_fdr = F, filtered_pairs_saved_in = "global_filtered_pairs"  ,
method = "spearman"
){
    if (exists(filtered_pairs_saved_in) == F) {
        global_filtered_pairs <<- data.frame()
        cat ("初始化了相应的变量,将剔除掉的变量对子保存在了 global_filtered_pairs 变量中. 这个名字是固定写死的,不要更改.")
    }

    part_1 <- left[, unique_id_colname] %>% as.character()
    part_2 <- right[, unique_id_colname] %>% as.character()
    part_3 <- confounding[, unique_id_colname] %>% as.character ()

    pooled_ID <- c(part_1, part_2, part_3) %>% unique()


    neworder <- pooled_ID %>% order()
    pooled_ID <- pooled_ID[neworder] %>% data.frame(. , stringsAsFactors = F)
    colnames(pooled_ID) <- unique_id_colname

    left <- merge (pooled_ID , left , by = c(unique_id_colname) , all.x = T)
    right <- merge (pooled_ID , right , by = c(unique_id_colname) , all.x = T)
    confounding <- merge (pooled_ID , confounding , by = c(unique_id_colname) , all.x = T)

    left <- left[, - 1]
    right <- right[, - 1]
    confounding <- confounding[, - 1]

    for (i in 1 : ncol (left)) {
        if (i == 1) { pooled_r_and_p <- data.frame()
            pooled_error <- data.frame ()
        }
        for (j in 1 : ncol (right)) {
            pair_1 <- colnames (left)[i]
            pair_2 <- colnames (right)[j]

            temp_frame <- cbind (left[, i], right[, j], confounding)
            temp_frame <- na.omit(temp_frame)
            temp_frame <- apply (temp_frame , 2, function (each_col){
                result <- each_col %>% as.character () %>% as.numeric()
            })
            IS_error <<- F
            if (T) {
                tryCatch(each_line <- ppcor::pcor.test(x = temp_frame[, 1], y = temp_frame[, 2],
                z = temp_frame[, 3 : ncol(temp_frame)], method = c(method)),
                error = function(e) {
                    IS_error <<- T
                    # cat (i , ' ', pair_1 , j, ' ', pair_2 , "的偏相关系数无法被计算\n")

                    each_filtered_pair <- data.frame (pair_1 = pair_1, pair_2 = pair_2)
                    global_filtered_pairs <<- rbind (global_filtered_pairs, each_filtered_pair)
                })
            }

            if (IS_error == F) {
                each_line <- data.frame (pair_1 = pair_1 , pair_2 = pair_2 , each_line , stringsAsFactors = F)
                colnames(each_line) <- c("pair_1", "pair_2", "r", "p", "statistic", "n", "gp", "Method")
                pooled_r_and_p <- rbind (pooled_r_and_p , each_line)
            }
        }
    }

    IS_filtered <- is.na (pooled_r_and_p[, "p"]) == T
    temp_filtered_data <- pooled_r_and_p[IS_filtered  , c("pair_1", "pair_2")]
    global_filtered_pairs <<- rbind (global_filtered_pairs , temp_filtered_data)
    pooled_r_and_p <- pooled_r_and_p[! IS_filtered,]

    if (mannual_n_for_fdr > 0) {n <- mannual_n_for_fdr} else {
        n <- nrow(pooled_r_and_p)}
    return (pooled_r_and_p)
}

configData <- read.csv(opt$config, header = F, stringsAsFactors = F) %>%
    set_colnames(c("arg", "value")) %>%
    column_to_rownames("arg")

corP <- configData["cor.p", "value"] %>%
as.numeric()
corFdr <- configData["cor.fdr", "value"] %>%
as.numeric()
cor <- configData["cor.coe", "value"] %>%
as.numeric()

sampleIds <- read.csv(opt$g, header = T, stringsAsFactors = F) %>%
    select(c("SampleID", "ClassNote")) %>%
    .$SampleID

superClassData <- read.csv(opt$i, header = T, stringsAsFactors = F) %>%
select(c("Metabolite", "Class"))

diffData <- read_csv("../../potential/Markers.csv")

if (nrow(diffData) == 0) {
    quit(status = 0)
}

diffNames <- diffData %>%
.$Metabolite

orignalData <- read.csv(opt$i, header = T, stringsAsFactors = F, check.names = F) %>%
    select(- c("HMDB", "KEGG", "Class")) %>%
    filter(Metabolite %in% diffNames) %>%
    as.data.frame()

sampleInfo <- read.csv(opt$g, header = T, stringsAsFactors = F) %>%
select(c("SampleID", "ClassNote"))
classNotes <- sampleInfo %>%
.$ClassNote %>%
unique()
groupNames <- c(classNotes, paste0(classNotes, collapse = "_vs_"))

for (groupName in groupNames) {
    groups <- unlist(strsplit(groupName, "_vs_"))
    sampleIds <- sampleInfo %>%
        filter(ClassNote %in% groups) %>%
        .$SampleID

    data <- orignalData %>%
        select(c("Metabolite", sampleIds)) %>%
        gather("ID", "Value", - Metabolite) %>%
        spread(Metabolite, "Value")

    parent <- str_replace_all(groupName, "_vs_", "_and_")
    createWhenNoExist(parent)
    files <- list.files(opt$et, full.names = T)
    for (file in files) {
        dirName <- basename(file) %>%
        file_path_sans_ext()
        finalParent <- paste0(parent, "/", dirName)

        extraData <- read_csv(file) %>%
            rename(Metabolite = 1) %>%
            gather("ID", "Value", - Metabolite) %>%
            spread(Metabolite, "Value")

        confounderData <- read_csv(opt$confounder) %>%
            rename(Metabolite = 1) %>%
            gather("ID", "Value", - Metabolite) %>%
            spread(Metabolite, "Value")

        dataIds <- data$ID
        extraIds <- extraData$ID
        confounderIds<-confounderData$ID
        interIds <- Reduce(intersect, list(dataIds, extraIds, confounderIds))

        extraData <- extraData %>%
            filter(ID %in% interIds) %>%
            arrange(factor(ID, levels = interIds))

        confounderData <- confounderData %>%
            filter(ID %in% interIds) %>%
            arrange(factor(ID, levels = interIds))

        inData<-data %>%
            filter(ID %in% interIds) %>%
            arrange(factor(ID, levels = interIds))

        print("=inData=")
        print(head(inData))
        if(nrow(inData)<2){
            next
        }

        pcorData <- zd_pcor(left = inData, right = extraData, confounding = confounderData)

        allData <- tibble(Node1 = pcorData$pair_1, Node2 = pcorData$pair_2, r = pcorData$r, P = pcorData$p) %>%
            mutate_at(vars(c("r")), function(x){
                ifelse(is.na(x), 0, x)
            }) %>%
            mutate_at(vars(c("P")), function(x){
                ifelse(is.na(x), 1, x)
            }) %>%
            mutate(FDR = p.adjust(P, method = "fdr"))

        if(nrow(allData)==0){
            next
        }

        createWhenNoExist(finalParent)
        corData <- allData %>%
            select(c("Node1", "Node2", "r")) %>%
            spread(Node1, "r") %>%
            mutate_at(vars(- c("Node2")), function(x){
                ifelse(is.na(x), 0, x)
            }) %>%
            rename(` ` = Node2)
        write.csv(corData, paste0(finalParent, "/Inter_r_Matrix.csv"), quote = T, row.names = F)

        pData <- allData %>%
            select(c("Node1", "Node2", "P")) %>%
            spread(Node1, "P") %>%
            mutate_at(vars(- c("Node2")), function(x){
                ifelse(is.na(x), 1, x)
            }) %>%
            rename(` ` = Node2)
        write_csv(pData, paste0(finalParent, "/Inter_P_Matrix.csv"))

        fdrData <- allData %>%
            select(c("Node1", "Node2", "FDR")) %>%
            spread(Node1, "FDR") %>%
            mutate_at(vars(- c("Node2")), function(x){
                ifelse(is.na(x), 1, x)
            }) %>%
            rename(` ` = Node2)
        write.csv(fdrData, paste0(finalParent, "/Inter_FDR_Matrix.csv"), quote = T, row.names = F)

        edgeData <- allData %>%
            filter(P < corP & FDR <= corFdr & abs(r) > cor) %>%
            filter(Node1 != Node2) %>%
            mutate(distName = {
                Node1 %>%
                map2_chr(Node2, function(x, y){
                    vec <- c(x, y) %>%
                    sort()
                    str_c(vec, collapse = ";")
                })
            }) %>%
            distinct(distName, .keep_all = T) %>%
            select(- c("distName"))%>%
            arrange(desc(abs(r)))
        write_csv(edgeData, paste0(finalParent, "/Network_Edges_for_Cytoscape.csv"))

        nodes <- unique(c(edgeData$Node1, edgeData$Node2))
        infoData <- tibble(Node = nodes) %>%
            mutate(Size = {
                Node %>%
                map_int(function(x){
                    edgeData %>%
                        filter(Node1 == x | Node2 == x) %>%
                        nrow()
                })
            }) %>%
            left_join(superClassData, by = c("Node" = "Metabolite")) %>%
            rename(Type = Class) %>%
            mutate_at(vars("Type"), function(x){
                replace_na(x, "Others")
            })
        write_csv(infoData, paste0(finalParent, "/Network_Nodes_for_Cytoscape.csv"))
    }
}




